function [dresults, IVDtilda] = estimate_rcnl( ...
    Pi0,           ...
    rho0,          ...
    s_jt,          ...
    p_jt,          ...
    prodid,       ...
    xlinear,       ...
    xnonlin,       ...
    zdemand,       ...
    zinside,       ...
    fe,            ...
    market,        ...
    dfull,         ...
    cluster,       ...
    model,         ...
    par_rest,      ...
    coefs_sim,     ...
    tolerance,     ...
    varargin       ...
)

    p = inputParser;
    addParameter(p, 'rhoTransform', 1);
    addParameter(p, 'rhoUnbounded', 0);
    addParameter(p, 'linearGMM',    0);
    parse(p, varargin{:});

    tic;
    options = optimoptions('fmincon',     ...
        'Display', 'iter-detailed',       ...
        'MaxIterations',150,              ...
        'Algorithm', 'interior-point',    ...
        'StepTolerance', 1e-10,           ...
        'SpecifyObjectiveGradient', true, ...
        'OptimalityTolerance',1e-8,       ...
        'CheckGradients',false            ...
    );


    rho_lb   = -Inf;
    rho_ub   = 8;
    rho_zero = -realmax('double');
    rho_transform = 0.95;
    rho0     = -log(rho_transform/rho0-1);
    if rho0  == -Inf
        rho0 = -realmax('double');
    end
    
    lb  = [zeros(1, length(Pi0)) rho_zero];
    ub  = [zeros(1, length(Pi0)) rho_zero];
    
    lb1 = [-Inf(1, length(Pi0)) rho_lb];
    ub1 = [Inf(1, length(Pi0)) rho_ub];
    
    derdefined = par_rest;
    
    if strcmp(model,'RCNL')
        derdefined  = [derdefined,4];
        lb(end)     = rho_lb;
        ub(end)     = rho_ub;
    end
    
    lb(par_rest) = lb1(par_rest);
    ub(par_rest) = ub1(par_rest);
    
    if strcmp(model,'RCNL_f')
        lb  = [Pi0' rho0];
        ub  = [Pi0' rho0];
    end    

    % -----------------------------------------
    % Initialize mean value for share inversion
    % -----------------------------------------

    cdid    = grp2idx(market);
    cdindex = splitapply(@(s) s(end), (1:size(s_jt, 1))', cdid);
    inshr   = accumarray(cdid, s_jt);
    inshr   = inshr(cdid);
    outshr  = 1.0 - inshr;
    logodds = log(s_jt) - log(outshr);

    expmval_mat = [tempname(), '.mat'];
    expmval0 = exp(logodds);
    save(expmval_mat, 'expmval0');

    IVDtilda    = zdemand - fe * ((fe' * fe) \ (fe' * zdemand));
    nonzerocols = (max(abs(IVDtilda))>0.00001); % remove columns with zeros
    IVDtilda    = IVDtilda(:,nonzerocols);
    invADtilda1 = inv(IVDtilda' * IVDtilda/length(cluster));
    dparams     = [Pi0; rho0];
    
    IVDtilda_alt      = zinside - fe * ((fe' * fe) \ (fe' * zinside));
    nonzerocols       = (max(abs(IVDtilda_alt))>0.00001); % remove columns with zeros
    IVDtilda_alt      = IVDtilda_alt(:,nonzerocols);
    invADtilda1_alt   = inv(IVDtilda_alt' * IVDtilda_alt/length(cluster));

    modelnums = '';
    for i = par_rest
        modelnums = [modelnums,int2str(i)];
    end

    disp(' ')
    fprintf("Estimating model %s%s\n", model, modelnums);

    % ----------
    % First step
    % ----------

    if ~strcmp(model,'RCNL_f')

    f_obj = @(dparams)rcnl_gmmobj_demand( ...
        dparams,           ...
        s_jt,              ...
        xlinear,           ...
        xnonlin,           ...
        cdindex,           ...
        dfull,             ...
        IVDtilda,          ...
        invADtilda1,       ...
        IVDtilda_alt,      ...
        invADtilda1_alt,   ...
        fe,                ...
        expmval_mat,       ...
        1,                 ...
        rho_transform,     ...
        tolerance          ...
    );
    f_obj(dparams);

    [dparams_1, ~, ~, ~] = fmincon(f_obj, dparams, [], [], [] ,[], lb, ub, [], options);

    [~, ~, ~, xi_1] = f_obj(dparams_1);

    % Optimal cluster re-weighting
    [S_1, ~, ~, W_new_1] = rcnl_dvarcov( ...
        dparams_1,         ...
        s_jt,              ...
        xlinear,           ...
        xnonlin,           ...
        cdindex,           ...
        dfull,             ...
        IVDtilda,          ...
        invADtilda1,       ...
        IVDtilda_alt,      ...
        invADtilda1_alt,   ...
        fe,                ...
        xi_1,              ...
        cluster,           ...
        expmval_mat,       ...
        rho_transform,     ...
        derdefined,        ...
        tolerance          ...
    );
    invADtilda2     = W_new_1; 
    invADtilda2_alt = invADtilda1_alt;

    if (strcmp(model,'RCNL') && length(par_rest)==3) 
        ncols = size(IVDtilda,2);
        invADtilda2     = inv(S_1(1:ncols,1:ncols));
        invADtilda2_alt = invADtilda2;
    end
    
    dhours_1 = toc / 60 / 60;
    disp(['Elapsed (total): ' num2str(dhours_1)]);


    % -----------
    % Second step
    % -----------

    % Re-estimate using optimal weights
    f_obj = @(dparams)rcnl_gmmobj_demand( ...
        dparams,           ...
        s_jt,              ...
        xlinear,           ...
        xnonlin,           ...
        cdindex,           ...
        dfull,             ...
        IVDtilda,          ...
        invADtilda2,       ...
        IVDtilda_alt,      ...
        invADtilda2_alt,   ...
        fe,                ...
        expmval_mat,       ...
        1,                 ...
        rho_transform,     ...
        tolerance          ...
    );

    [dparams_2, ~, ~, ~] = fmincon(f_obj, dparams_1, [], [], [] ,[], lb, ub, [], options);

    dhours_2 = toc / 60 / 60;
    [~, ~, theta_2, xi_2] = f_obj(dparams_2);
    disp(['Elapsed (total): ' num2str(dhours_2)]);

    else % fix alpha for estimation with fixed parameters

    dparams_2 = dparams;

    Pi_2  = dparams(1:end-1);
    rho_2 = dparams(end);
    if rho_transform
        rho_2 = rho_transform * (1/(1+exp(-rho_2)));
    end

    [delta_2, ~] = rcnl_meanval(Pi_2, rho_2, s_jt, xnonlin, cdindex, dfull, expmval_mat, 1, tolerance);
    yfe          = delta_2 - xlinear(:,1) * coefs_sim.alpha;
    theta_2      = inv(fe' * fe) * fe' * yfe;
    theta_2      = [coefs_sim.alpha;theta_2]; 
    xi_2         = delta_2 - xlinear*theta_2;


    end

    delta = NaN;
    
    % Estimate elasticity

    [ownElas, derMat]   = rcnl_der( ...
        theta_2(1),         ...
        dparams_2(1:end-1), ...
        dparams_2(end),     ...
        s_jt,               ...
        p_jt,               ...
        xnonlin,            ...
        cdindex,            ...
        dfull,              ...
        market,             ...
        prodid,             ...
        delta,              ...
        expmval_mat,        ...
        0,                  ...
        rho_transform,      ...
        tolerance           ...
    );


    % --------------
    % Gather results
    % --------------

    dresults           = struct;
    dresults.ownElas   = ownElas;
    dresults.Pi        = dparams_2(1:end-1);
    dresults.theta     = theta_2;
    dresults.rho       = dparams_2(end);
    dresults.xi        = xi_2;
    dresults.derMat    = derMat;
    if length(par_rest)<3
        dresults.derMat                    = NaN;
    end

    if rho_transform
        dresults.rho              = rho_transform * (1/(1+exp(-dparams_2(end))));
    end

end
